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General ized-a  Time  Integration  Solutions 
for  Hanging  Chain  Dynamics 


Jason  I.  Gobat1;  Mark  A.  Grosenbaugh2;  and  Michael  S.  Triantafyllou3 


Abstract:  In  this  paper,  we  study  numerically  the  two-  and  three-dimensional  nonlinear  dynamic  response  of  a  chain  hanging  under  its 
own  weight.  Previous  authors  have  employed  the  box  method,  a  finite-difference  scheme  popular  in  cable  dynamics  problems,  for  this 
purpose.  The  box  method  has  significant  stability  problems,  however,  and  thus  is  not  well  suited  to  this  highly  nonlinear  problem.  We 
illustrate  these  stability  problems  and  propose  a  new  time  integration  procedure  based  on  the  generalized-a  method.  The  new  method 
exhibits  superior  stability  properties  compared  to  the  box  method  and  other  algorithms  such  as  backward  differences  and  trapezoidal  rule. 
Of  four  time  integration  methods  tested,  the  generalized-a  algorithm  was  the  only  method  that  produced  a  stable  solution  for  the 
three-dimensional  whirling  motions  of  a  hanging  chain  driven  by  harmonic  linear  horizontal  motion  at  the  top. 

DOI:  10.1061/(ASCE)0733-9399(2002)  128:6(677) 
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Introduction 

The  dynamics  of  a  chain  hanging  under  its  own  weight  is  a  classic 
problem  in  mechanics.  Two  of  the  more  interesting  aspects  of  the 
problem  are  the  simultaneous  presence  of  both  high-  and  low- 
tension  regimes  in  the  chain  and  the  unstable  nature  of  large  am¬ 
plitude  motions.  Triantafyllou  and  Howell  (1993)  and  Howell  and 
Triantafyllou  (1993)  considered  both  of  these  phenomena  using  a 
combination  of  analytic,  numerical,  and  experimental  results. 
They  observed  that  the  stability  of  the  response  in  a  harmonically 
driven  system  is  strongly  dependent  on  the  frequency  and  ampli¬ 
tude  of  the  excitation. 

The  numerical  model  that  they  employed  was  based  on  a 
finite-difference  scheme  known  as  the  box  method.  This  method 
was  first  applied  to  a  cable  dynamics  problem  by  Ablow  and 
Schechter  (1983).  Because  the  box  method  is  an  implicit  scheme, 
box  method  solutions  for  the  classical  cable  dynamics  equations 
are  singular  when  the  tension  goes  to  zero  anywhere  on  the  cable. 
Howell  and  Triantafyllou  (1993)  removed  this  singularity  by  add¬ 
ing  bending  stiffness  to  the  governing  equations,  thus  providing  a 
mechanism  to  propagate  energy  in  the  presence  of  zero  tension 
(Burgess  1993).  For  small  values  of  artificial  bending  stiffness 
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this  modification  stabilized  the  numerical  solution  with  no  loss  of 
accuracy  compared  to  experimental  results. 

The  box  method  is  popular  because  it  is  second-order  accurate 
in  both  space  and  time  and  is  relatively  easy  to  implement.  Be¬ 
cause  the  box  method  preserves  the  frequency  content  of  the  so¬ 
lution  across  all  frequencies,  however,  it  has  the  disadvantage  of 
relatively  poor  stability  in  its  temporal  discretization.  In  a  nonlin¬ 
ear  problem,  spurious  high-frequency  content  can  cause  numeri¬ 
cal  instabilities,  and  thus,  it  is  desirous  that  a  temporal  integration 
scheme  should  be  numerically  dissipative  at  high  frequencies. 
Koh  et  al.  (1999)  addressed  this  shortcoming  of  the  box  method 
by  replacing  the  box  method’s  temporal  integration  scheme  with 
backward  differences.  They  preserved  the  box  method’s  straight¬ 
forward  and  easy  to  implement  spatial  discretization.  Backward 
differences  have  also  been  used  by  Chatjigeorgiou  and  Mavrakos 
(1999)  and  Chiou  and  Leonard  (1991)  in  conjunction  with  spatial 
discretizations  based  on  collocation  and  direct  integration,  respec¬ 
tively.  The  scheme  is  only  first-order  accurate,  but  is  very  stable 
because  it  has  strong  numerical  dissipation  at  high  frequencies. 

Another  temporal  integration  scheme  that  has  been  used  in 
cable  dynamics  applications  is  the  generalized  trapezoidal  rule 
(Sun  et  al.  1994).  This  scheme  offers  controllable  numerical  dis¬ 
sipation,  but  is  second-order  accurate  only  in  its  least  dissipative 
form.  Thomas  (1993)  compared  three  historically  popular  algo¬ 
rithms  from  the  structural  dynamics  community,  Newmark, 
Houbolt,  and  Wilson-0,  for  use  in  mooring  dynamics  problems. 
His  conclusion  was  that  Houbolt  was  the  best  choice  of  the  three. 
Other  authors,  however,  have  noted  that  Houbolt  has  an  undesir¬ 
able  amount  of  low-frequency  dissipation  (Chung  and  Hulbert 
1994;  Hughes  1987). 

Turning  to  the  more  recent  structural  dynamics  literature, 
Gobat  and  Grosenbaugh  (2001)  proposed  replacing  the  box  meth¬ 
od’s  temporal  integration  with  the  generalized-a  method  devel¬ 
oped  for  the  second-order  structural  dynamics  problem  by  Chung 
and  Hulbert  (1993).  This  algorithm  has  the  advantages  of  control¬ 
lable  numerical  dissipation,  second-order  accuracy,  and  straight¬ 
forward  adaptation  to  the  first-order  nonlinear  cable  dynamics 
problem.  Through  appropriate  choices  of  parameters,  the  method 
can  also  reproduce  the  spectral  properties  of  several  other  algo- 
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rithms  including  the  box  method,  backward  differences,  and  trap¬ 
ezoidal  rule.  This  latter  property  makes  it  a  particularly  conve¬ 
nient  choice  for  the  type  of  comparative  study  undertaken  herein. 

The  analyses  of  the  box  and  generalized-a  methods  from 
Gobat  and  Grosenbaugh  (2001)  are  summarized  below.  The  per¬ 
formance  of  the  new  algorithm  is  studied  by  comparison  to  ana¬ 
lytic  and  experimental  results  for  the  free  and  forced  response  of 
the  hanging  chain.  Throughout  the  analyses,  comparisons  are  also 
made  to  trapezoidal  rule  and  backward  difference  solutions. 


Analysis  of  Box  Method 

The  governing  equations  for  a  cable  or  chain  can  be  written  as  a 
system  of  partial  differential  equations  of  the  form  (Howell  1992) 

dY  dY 

M  — +  K— -f  F(Y,s,f)  =  0  (1) 

where  Y=  vector  of  ^-dependent  variables,  M  and  K 
=  coefficient  matrices,  and  F = force  vector.  The  independent 
variables  are  s,  the  Lagrangian  coordinate  measuring  length  along 
the  unstretched  cable,  and  f,  time.  Howell  and  Triantafyllou 
(1993)  used  the  box  method  to  discretize  Eq.  (1).  In  the  box 
method  the  discrete  equations  are  written  using  what  look  like 
traditional  backward  differences  in  both  space  and  time,  but  be¬ 
cause  the  discretization  is  applied  on  the  half-grid  points  with 
spatial  and  temporal  averaging  of  adjacent  grid  points,  the  method 
is  second-order  accurate.  The  result  is  a  four-point  average  cen¬ 
tered  around  the  half-grid  point. 

The  stability  of  the  box  method  can  be  analyzed  by  consider¬ 
ing  an  equivalent  linear,  single  degree-of-freedom  system  in  se¬ 
midiscrete  form.  This  approach  separates  the  spatial  and  temporal 
discretizations  into  distinct  procedures.  For  each  of  the  n  —  1  spa¬ 
tial  half-grid  points  between  the  n  nodes  a  set  of  N  discrete  equa¬ 
tions  is  assembled.  Combining  these  N(n  —  1 )  equations  with  N 
equations  describing  the  boundary  conditions  yields  the  semidis¬ 
crete  equation  of  motion  for  all  of  the  dependent  variables  at  all 
of  the  nodes  as  (Gobat  and  Grosenbaugh  2001) 

MY+KY+F=0  (2) 

The  tilde  over  the  matrices  signifies  that  these  are  now  dis¬ 
cretized,  assembled  quantities.  The  single  degree-of-freedom,  lin¬ 
ear,  homogeneous  analog  of  Eq.  (2)  is 

y  +  (oy  =  0  (3) 

Applying  the  box  method’s  temporal  discretization  to  Eq.  (3) 


yields 

y*+yi~1  +  <o(yz+y'  l)=0 

(4) 

where 

(5) 

Rearranging  Eq.  (5)  gives  the  recursion  relationships 

(6) 

.  Af  .  ,  .  , 

y,=Y^y,+y‘  ')+y'  1 

(7) 

Substituting  each  of  the  recursion  relationships  separately  into  Eq. 
(4),  we  can  write  equations  for  yl  and  yl  in  matrix  form  as 


/ 

y 


2-o)A  t 
2  +  coA  t 


L  2  +  coAf 


(8) 


The  2X2  matrix  on  the  right-hand  side  of  Eq.  (8)  is  the  am¬ 
plification  matrix.  Spectral  radius  p  of  this  matrix,  defined  as 

p=max(|XI|,|\2|)  (9) 

governs  the  growth  or  decay  of  the  solution  from  one  time  step  to 
the  next  (Hughes  1987).  X12=  eigenvalues  of  the  amplification 
matrix.  For  1,  the  solution  will  remain  steady  or  decay  and  is 
said  to  be  stable.  For  p>  1,  the  solution  will  grow  and  is  said  to 
be  unstable.  For  the  box  method, 


2  —  G)Ar 
Xl  =  2  +  wAr 


(10) 


X2=  —  1  (11) 

and  the  spectral  radius  is  unity  (and  the  scheme  is  stable)  for  all 
values  of  to  and  At. 

In  spite  of  this  unconditional  stability,  however,  the  box 
method  has  three  significant  problems.  The  first  problem  is  illus¬ 
trated  by  considering  the  update  equation  for  yl  written  in  the 
form 


y« 


2  —  toAr  \ 

2  +  o)A  t)y 


(12) 


As  to  At  goes  to  infinity  this  becomes 


yi=-yi  1 


(13) 


This  is  the  phenomenon  known  as  Crank-Nicholson  noise, 
whereby  the  high-frequency  components  of  the  solution  oscillate 
with  every  time  step.  A  second,  related,  problem  is  that  the  spec¬ 
tral  radius  is  constant  at  unity.  An  artifact  of  the  spatial  discreti¬ 
zation  process  is  that  at  some  point  the  high-frequency  (or  equiva¬ 
lently,  high-spatial  wave-number)  components  of  the  solution  are 
not  well  resolved  and  the  numerical  solution  is  inaccurate.  For 
this  reason  it  is  desirous  to  have  numerical  dissipation  in  a 
scheme  such  that  the  spectral  radius  is  less  than  unity  for  increas¬ 
ing  values  of  to  A  t.  The  box  method  has  no  numerical  dissipation. 
Finally,  Hughes  (1977)  cites  a  problem  with  averaging  schemes  in 
general  as  applied  to  nonlinear  problems.  For  the  nonlinear  single 
degree-of-freedom  case,  Eq.  (4)  can  be  written  as 

y'+y'"1  +  <o*y2  +  a)I-1y1-1  =  0  (14) 


The  update  equation  for  y\  Eq.  (12),  becomes 


/  = 


2-co‘”1Ar\ 

2  +  w‘Ar  K 


(15) 


and  the  stability  becomes  conditional  as  parameter  w  changes 
with  time.  The  practice  suggested  by  Hughes  (1977)  for  avoiding 
this  problem  is  to  use  an  averaged  value  of  o>,  i.e., 

y'  +  y'"1-^ - ^ - )  (y' 1)  =  0  (16) 


Generalized-a  Method 

Given  the  stability  problems  associated  with  the  box  method, 
Gobat  and  Grosenbaugh  (2001)  proposed  replacing  the  temporal 
integration  with  Chung  and  Hulbert’s  (1993)  generalized-a 
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Table  1 .  Algorithms  Included  in  Generalized-a  Method 


Algorithm 

7 

a* 

1st  order  problem 

2nd  order  problem 

Box  method 

\ 

2 

1 

2 

Ablow  and  Schechter  (1983) 

Backward  differences 

1 

0 

0 

Koh  et  al.  (1999) 

Generalized  trapezoidal 

[iU 

0 

0 

Sun  et  al.  (1994) 

Newmark  (1959) 

Cornwell  and  Malkus 

i 

j- a 

a 

0 

Cornwell  and  Malkus  (1992) 

Hilber  et  al.  (1977) 

WBZ-a 

0 

a 

Wood  et  al.  (1981) 

method.  The  generalized-a  method  is  a  reasonably  complete  fam¬ 
ily  of  algorithms  that  is  second-order  accurate,  has  controllable 
numerical  dissipation,  and  offers  a  clear  approach  to  coefficient 
averaging  for  the  nonlinear  problem.  Following  Chung  and  Hul- 
bert’s  development  of  the  generalized-a  method  for  second-order 
equations,  semidiscrete  Eq.  (2)  becomes 

( 1  -  aJMY1  +  amMY/“ 1  +  ( 1  -  a*)KYl  +  a*KY‘ ~ 1 

+  (l-a*)F‘  +  a*F'~1  =  0  (17) 

The  difference  equation  is  the  same  as  for  the  generalized  trap¬ 
ezoidal  rule  (Hughes  1987), 

Y«  =  Y‘- 1  +  A/[(  1  -  y )  Y,_  1  +  7  Y1']  (18) 

The  three  parameter  family  of  algorithms  given  by  Eqs.  (17)  and 
(18)  defines  the  generalized-a  method  for  the  first-order  semidis¬ 
crete  problem.  The  method  is  second-order  accurate  if 

am~ak  +  'y~2  (19) 

From  the  eigenvalues  of  the  amplification  matrix,  the  stability 
requirement  is 

(20) 


Requiring  second-order  accuracy  according  to  Eq.  (19)  and  forc¬ 
ing  the  eigenvalues  of  the  amplification  matrix  to  be  equal  as 
o)A/— ►«>  to  prevent  bifurcation,  yields  formulas  for  <xk  and  am  as 
a  function  of  X00  only 

X06  3X°°  + 1 

^  2\°°  —  2 

This  yields  a  second-order  accurate  algorithm  in  which  the  only 
parameter  is  the  eigenvalue  (or  spectral  radius)  at  infinity. 

Algorithms  that  can  be  obtained  through  various  choices  of 
a* ,  am ,  7,  and  X00  are  listed  in  Table  1.  Spectral  radii  of  some  of 
these  algorithms  are  shown  in  Fig.  1.  Note  that  taking  X00 
e  [0,1)  as  the  basis  for  the  spectral  radius  results  in  a  different  set 
of  algorithms  than  X*  e  [ -  1,0].  For  p°°=  1  the  only  option  is  the 
negative  eigenvalue  and  this  results  in  the  box  method.  A  nondis- 
sipative  algorithm  with  X“=  + 1  cannot  be  achieved. 

In  applying  the  generalized-a  method  to  the  nonlinear  problem 
we  must  choose  the  time  point  at  which  we  will  evaluate  M,  K, 
and  F.  A  natural  choice,  consistent  with  the  practice  suggested  by 
Hughes  (1977)  for  nonlinear  first-order  problems  and  exemplified 
by  Eq.  (16),  is  provided  by  the  temporal  averaging  of  terms  that  is 
already  a  part  of  the  method.  At  time  step  i  Eq.  (17)  becomes 


Fig.  1 .  Spectral  radii  of  generalized-a  family  algorithms 
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( 1  -  am)M‘‘-a"Y‘  +  amMI'-“»Yi  - 1  +  ( 1  -  a^KT^Y'' 

+  oikK-^Y- 1  +  ( 1  -  afc)F  +  ahP~ 1  =  0  (22) 

where  the  averaged  coefficient  matrices  are  defined  as 

Mi“a«  =  (l  -am)Mf  +  a JMT 1  (23) 

K*-«t=(  i  -  a*)R*  +  1  (24) 

This  scheme  has  been  implemented  in  a  computer  program  for 
two-  and  three-dimensional  simulations  of  cable  dynamics  (Gobat 
and  Grosenbaugh  2000).  At  each  time  step,  Eq.  (22)  is  solved 
using  a  Newton-Raphson  procedure.  The  solution  from  the  pre¬ 
vious  time  step  (or  the  static  solution  at  the  initial  time  step) 
serves  as  the  initial  guess  in  the  nonlinear  iterations.  Because  of 
this,  the  ultimate  success  of  the  solution  is  dependent  on  both  the 
stability  of  the  time  integration  and  on  the  ability  of  the  nonlinear 
solver  to  converge  on  a  solution  at  time  step  i  given  an  initial 
guess  based  on  the  solution  at  time  step  i  —  1 .  To  improve  conver¬ 
gence  the  program  implements  an  adaptive  time  stepping  scheme 
whereby  the  time  step  (the  distance  between  the  guess  at  i  —  1  and 
the  solution  at  i)  is  reduced  by  factors  of  10  at  any  spots  where 
the  solver  is  not  successful.  A  practical  limit  of  four  orders  of 
magnitude  below  the  base-line  time  step  is  set  to  prevent  the 
solution  from  proceeding  in  the  face  of  a  physical  or  numerical 
instability  unrelated  to  the  nonlinear  solution  procedure  (e.g., 
Crank-Nicholson  noise). 

All  of  the  numerical  solutions  that  follow  were  obtained  using 
this  program.  Thus,  the  box  method,  trapezoidal  rule,  and  back¬ 
ward  difference  results,  while  spectrally  equivalent  to  previous 
implementations,  may  be  more  stable  than  previous  solutions  be¬ 
cause  of  the  coefficient  averaging  scheme  in  Eq.  (22).  For  clarity, 
spectrally  equivalent  historical  names  are  retained  in  discussions 
of  comparative  algorithm  performance  that  follow. 


Application  to  Hanging  Chain  Problem 

The  performance  of  the  different  algorithms  that  can  be  imple¬ 
mented  with  the  generalized-a  family  is  studied  by  considering 
the  free  and  forced  response  of  the  hanging  chain  shown  in  Fig.  2. 
In  the  free-response  problem,  we  apply  a  small  initial  displace¬ 
ment  to  the  chain  and  then  at  time  t  =  0,  release  it.  The  dynamic 
response  of  the  chain  for  *>0  can  be  calculated  analytically  for 
the  small  motions  that  result.  In  the  forced  response  problem  we 
impose  a  sinusoidally  varying  horizontal  displacement  to  the  top 
of  the  chain  and  analyze  the  forced  response.  This  latter  problem 
was  studied  both  numerically  and  experimentally  by  Howell  and 
Triantafyllou  (1993). 


Free  Response  to  Initial  Displacement 

For  small  motions  and  an  inextensible  chain,  the  equation  of  mo¬ 
tion  is 


m 


d 

ds 


mgs 


dq 

ds 


(25) 


where  m  =  mass  per  length  of  the  chain,  q  =  transverse  displace¬ 
ment  of  the  chain,  g  =  acceleration  due  to  gravity,  and  5 
=  independent  coordinate  along  the  chain  with  5  =  0  at  the  free 
end.  Assuming  a  harmonic  solution  of  the  form 

q(syt)=q{s)[A  coscof+B  sin  cor]  (26) 

the  mode  shapes;  q(s ),  are  (Triantafyllou  et  al.  1986) 
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^Q(t) 


Fig.  2.  Definitions  for  hanging  chains  problems 


q(s)  =  c1J0[2uyJ^ 


+c2yq\  2 


(27) 


where  J0  and  Y0=  zero-order  Bessel  functions  of  the  first  and 
second  kind,  respectively.  The  requirement  that  the  solution  be 
finite  at  5  =  0  leads  to  the  elimination  of  the  Y0  term  and  the 
requirement  that  q(L)  =  0  leads  to  the  natural  frequencies,  o>. 
They  are  given  by  the  roots  of 


J  n  2(0 


j).. 


(28) 


The  complete  response  is  given  as  the  sum  of  the  response  in  all 
modes: 


q(,s,t)=Yi  Jo  2w« 

n=  1  \ 


[An  coscoi  +  Bn  sinu>r]  (29) 


The  coefficients  An  and  B„  are  determined  from  the  initial 
displacement,  q0(s ),  and  velocity,  qo(s)-  Given  4o(J) =  we  can 
immediately  determine  that  B„  =  0.  To  determine  A„  we  first  write 


q(s 


,0)=S  A„J0\2<*n^=qo(s) 


(30) 


Multiplying  both  sides  by  JQ(2vn\[s7g ),  integrating  from  5  -  0  to 
5  =  Ly  and  making  use  of  the  fact  that 


J  7°|  2o)n  2 wm  ^"^ds  — 


/o[2o)a 

yields  the  following  equation  for  An  : 

I  ft 

q o(-s)A)  2wn  V”j 

Jo  \  y  g 


0  for  (31) 


ds 


A  =- 


(32) 


ds 
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Fig.  4.  Power  spectra  peaks  of  response  of  free  end  of  chain  for  analytic  solution  and  for  six  variants  of  generalized-a  method  with  A  t 
=  0.001  s  and  n  =  50 
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non-dimensional  frequency 

Fig.  5.  Power  spectra  peaks  of  response  of  free  end  of  chain  for  analytic  solution  and  for  \x=  -  5,  At =0.001  s,  with  n  — 50,  200,  400,  and  800 
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Fig.  6.  Snapshots  of  chain  configuration  near  time  of  expected  intersection  for  six  algorithms  that  can  be  obtained  from  within  generalized-a 
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Fig.  7.  Snapshots  of  chain  configuration  near  time  of  expected  intersection  for  box  method  with  different  spatial  and  temporal  discretizations 


The  analytic  solution  was  computed  for  a  chain  released  from 
an  initial  catenary  configuration.  For  simplicity  all  of  the  model 
parameters  (mass  per  length,  gravity,  length)  were  set  to  unity. 
The  horizontal  force  applied  at  s  =  0  to  create  the  initial  deflection 
was  set  to  0.001  N.  All  of  the  integrals  for  the  analytic  solution 
were  computed  using  the  trapezoidal  rule  with  10,000  panels.  A 
400  s  time  series  of  the  response  at  the  free  end  was  constructed 
using  the  first  20  modes  of  the  analytic  solution.  The  analytic 
result  was  sampled  at  20  Hz  to  adequately  capture  the  response  up 
to  mode  20.  (The  natural  frequency  for  mode  20  is  approxi¬ 


mately  5  Hz.) 

Analytic  solutions  were  compared  to  numerical  simulation  re¬ 
sults  for  a  chain  released  from  the  same  initial  configuration.  For 
simulation  results  El  was  set  to  10“ 6  Nm2  and  EA  to  109  N.  This 
setting  for  El  corresponds  to  the  value  of  El*-EHmgl}  used  in 
Howell’s  (1992)  comparison  of  experimental  and  simulation  re¬ 
sults  and  in  the  simulations  of  the  forced  hanging  chain  problem 
that  follow.  The  results  from  Howell  demonstrated  that  this  value 
is  sufficient  to  stabilize  the  numerical  solution  in  the  presence  of 
zero  tension,  but  is  small  enough  as  to  have  a  negligible  effect  on 
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Fig.  8.  Snapshots  of  chain  configuration  near  time  of  expected  intersection  for  trapezoidal  rule  with  different  spatial  and  temporal  discretizations 
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the  accuracy  of  the  simulation  result  (based  on  comparisons  with 
experiment).  The  model  results  were  insensitive  to  changes  in  this 
value  of  at  least  an  order  of  magnitude. 

Because  the  primary  distinction  among  the  various  algorithms 
contained  within  the  generalized-a  family  is  the  amount  of  nu¬ 
merical  dissipation,  all  results  are  compared  in  the  frequency  do¬ 
main.  For  each  400  s  time  series,  power  spectra  of  the  response  at 
the  free  end  were  computed  using  nonoverlapping  256  point  fast 
Fourier  transforms.  For  clarity,  only  the  peaks  of  the  spectra  are 
plotted.  This  prevents  clutter  and  allows  for  a  comparison  of  the 
spectral  roll  off  of  each  of  the  algorithms  compared  to  the  roll  off 
from  the  analytic  solution. 

Fig.  3  shows  a  comparison  between  the  analytic  solution  and 
numerical  solutions  for  six  different  parametrizations  of  the 
generalized-a  method.  At  this  time  step,  At  =  0.01  s,  most  of  the 
algorithms  are  accurate  out  to  the  fifth  or  sixth  mode.  The  notable 
exception  is  the  first-order  accurate  backward  difference  solution, 
which  substantially  underestimates  the  response  even  in  the  first 
mode.  All  of  the  algorithms  show  a  marked  fall  off  from  the 
analytic  solution  at  higher  frequencies,  with  the  solutions  for  \x 
>0  showing  the  most  decay  and  the  trapezoidal  rule  appearing  to 
be  the  most  accurate. 

In  Fig,  1,  the  numerical  damping  of  most  algorithms  increases 
with  increasing  <*>At.  The  idea  that  we  should  see  less  numerical 
damping  at  a  fixed  frequency  with  a  decrease  in  At  is  illustrated 
in  Fig.  4,  which  shows  the  same  results  comparison  as  in  Fig.  3 
for  a  time  step  of  At= 0.001  s.  At  this  time  step  most  algorithms 
are  accurate  out  to  the  tenth  mode.  Only  backward  differences, 
which  due  to  its  first-order  accuracy  is  again  a  poor  solution  even 
at  very  low  frequencies,  and  A°°=0  are  worse  than  this. 

That  the  other  algorithms,  with  their  varying  levels  of  dissipa¬ 
tion,  have  converged  to  the  same  solution  suggests  that  the  re¬ 
maining  error  is  not  due  to  numerical  dissipation.  Fig.  5  shows  the 
comparison  for  four  cases  with  X°°  =  —  \  and  At— 0.001  s,  with  a 
varying  number  of  nodes.  As  the  node  density  is  increased,  the 
numerical  model  is  better  able  to  resolve  the  mode  shapes  asso¬ 
ciated  with  the  higher  frequencies.  At  n  =  800,  the  numerical  so¬ 
lution  is  in  agreement  with  the  analytic  solution  over  the  full 
range  of  the  analytically  computed  response. 

These  results  demonstrate  that  the  ability  of  the  model  to  ac¬ 
curately  resolve  high-frequency  response  is  dependent  on  tempo¬ 
ral  and  spatial  discretizations  and  on  the  numerical  dissipation  for 
a  given  algorithm.  Given  sufficient  temporal  and  spatial  resolu¬ 
tion,  most  of  the  algorithms  appear  ultimately  capable  of  accu¬ 
rately  calculating  the  free  response  of  the  swinging  chain.  Based 
on  its  better  accuracy  at  the  larger  time  step,  the  best  choice  of 
algorithm  for  this  problem  appears  to  be  the  trapezoidal  rule. 

Two-Dimensional  Forced  Response  to  imposed 
Motion 

The  forced  hanging  chain  problem  that  we  consider  was  studied 
by  Howell  and  Triantafyllou  (1993).  In  this  problem,  a  1.75-m- 
long  chain  is  suspended  from  an  actuator  which  imposes  a  sinu¬ 
soidally  varying  horizontal  linear  displacement,  Q(t),  to  the  top 
of  the  chain  (Fig.  2).  In  experiments,  Howell  and  Triantafyllou 
observed  that  the  free  end  of  the  chain  intersects  the  chain  above 
it  at  approximately  3.4  s. 

Fig.  6  shows  the  configuration  of  the  lower  portion  of  the 
chain  from  3.43  to  3.46  s  for  six  different  numerical  algorithms, 
all  with  n  =  100  and  Af  =  0.01  s.  The  box  method  and  trapezoidal 
rule  both  closely  match  the  experimental  result,  with  intersection 
occurring  by  the  3.43  s  mark.  For  the  other  algorithms  the  inter¬ 


section  occurs  later;  the  delay  in  the  time  of  intersection  is  pro¬ 
portional  to  the  amount  of  numerical  dissipation  in  the  algorithm. 
The  backward  differences  solution  is  again  the  worst;  the  chain 
never  intersects  itself.  Likewise  for  Ax  =  0,  though  it  comes  closer 
to  doing  so.  For  A°°  =  -0.7,  intersection  actually  happens  at  3.47 
s  and  for  A°°  =  -0.5,  at  3.50  s. 

The  situation  changes  somewhat  if  we  consider  the  effect  of 
temporal  and  spatial  discretization.  Fig.  7  shows  the  same  time 
points  for  versions  of  the  box  method  with  n=  100  or  200  and 
Af  =  0.01,  0.001,  and  0.0001  s.  In  this  case  we  see  that  increasing 
the  number  of  nodes  does  not  significantly  affect  the  solution, 
suggesting  that  n=  100  is  adequate  to  accurately  capture  the  re¬ 
sponse.  An  increase  in  temporal  resolution,  however,  from  Af 
=  0.01  s  to  Ar  =  0.001  s,  leads  to  a  delay  in  the  crossover  to  ap¬ 
proximately  3.46  s.  The  result  at  the  even  smaller  Af  =  0.0001  s 
confirms  that  the  solution  has  converged  at  these  smaller  time 
steps.  Fig.  8  shows  this  same  behavior  for  the  trapezoidal  rule. 
The  only  notable  difference  between  trapezoidal  rule  and  box 
method  solutions  is  the  better  smoothness  of  the  trapezoidal  rule 
solutions  at  Af=0.01  s. 

Similar  results  for  A°°=  -0.5  are  shown  in  Fig.  9.  In  this  case, 
the  solution  at  Af= 0.001  s  is  slightly  different  than  the  solutions 
from  the  trapezoidal  rule  and  the  box  method  at  the  3.46  s  point. 
The  solutions  for  Af= 0.0001  s  are  in  good  agreement  with  the 
converged  solutions  for  Af  =  0.001  s  in  Figs.  7  and  8.  A  notable 
difference  in  the  solutions  for  the  various  algorithms  does  appear 
between  3.5  and  4.0  s  (i.e.,  following  crossover).  Both  trapezoidal 
rule  and  box  method  solutions  required  significant  adaptation  of 
the  time  step  to  get  through  the  collapse  of  the  lower  portion  of 
the  chain  following  the  crossover.  The  enhanced  stability  of  solu¬ 
tions  with  A“=  —  0.5  allowed  for  a  smooth  numerical  solution  in 
this  region,  with  no  or  very  little  adaptation.  Without  experimen¬ 
tal  verification,  however,  we  cannot  say  if  the  X°°=  —0.5  solution 
is  accurate. 

At  sufficiently  small  time  steps  and  adequate  spatial  resolu¬ 
tion,  all  three  algorithms:  box  method,  trapezoidal  rule,  and  X00 
=  -0.5,  provide  accurate  solutions.  Trapezoidal  rule  is  the  best 
choice  in  terms  of  the  computational  cost  of  accuracy,  where  cost 
is  measured  in  terms  of  time  step.  As  indicated,  however,  in  re¬ 
gions  where  the  solution  becomes  numerically  unstable  some  nu¬ 
merical  dissipation  may  be  necessary  to  obtain  a  solution.  This 
suggests  a  trade  off  between  optimizing  the  time  step  for  accuracy 
and  optimizing  the  algorithm  for  stability. 

Three-Dimensional  Forced  Response  to  Imposed 
Motion 

In  order  to  further  explore  these  trade  offs,  three-dimensional 
simulations  were  conducted  to  explore  the  behavior  of  the  solu¬ 
tions  beyond  the  time  when  the  chain  crosses  over  itself.  Howell 
(1992)  noted  that  out-of-plane  motions  of  the  experimental  chain 
only  become  significant  after  this  point.  The  simulations  were 
conducted  with  a  small  initial  out-of-plane  force  applied  at  the 
free  end  to  promote  the  initiation  of  out-of-plane  motion.  This 
models  the  inevitable  presence  of  small  disturbing  forces  that 
produce  instabilities  in  the  two-dimensional  motion  and  eventu¬ 
ally  lead  to  a  fully  three-dimensional  (3D)  response. 

Table  2  lists  the  observed  time  of  the  chain  crossing  over  itself 
and  the  total  running  time  (out  of  a  possible  10  s  simulation)  of 
the  simulation  before  failure.  Only  solutions  for  —  0.4*sX°°^ 
—  0.2  ran  for  the  full  10  s  and  resulted  in  an  accurate  crossover 
prediction.  At  At =0.0 Is,  the  numerically  stable  solutions  (at 
X“=— 0.3  and  X°°=  —  0.2)  were  less  accurate  than  the  two- 
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Fig.  9.  Snapshots  of  chain  configuration  near  time  of  expected  intersection  for  \x  —  —  \  with  different  spatial  and  temporal  discretizations 


Table  2.  Crossover  Time  and  Total  Runtime  in  Three-Dimensional  Simulations 
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Fig.  10.  Out-of-plane  motion  of  free  end  of  hanging  chain  for  ~0.3  and  X  —  0.2  and  At  —  0.01  s  and  A t  0.001  s 
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dimensional  simulations  for  X 00  =-0.5  at  this  same  time  step. 
This  is  consistent  with  the  observation  that  as  damping  increases 
the  crossover  time  is  delayed.  Also  consistent  with  the  two- 
dimensional  results  is  the  convergence  to  a  prediction  of  3.46  s 
with  an  increase  in  temporal  resolution  to  At =0.001  s. 

The  numerical  stability  of  results  for  X"  =  0.0  and  X°°  =  0.1  at 
Ar=0.01  s,  but  not  at  At  =  0.001  s,  illustrates  the  dependence  of 
the  stability  on  the  frequency  content  of  the  response,  the  time 
step,  and  the  damping  properties  of  the  algorithm.  Because  the 
spectral  radii  in  Fig.  1  all  initially  decrease  with  the  product  go  At, 
a  decrease  in  A t  at  a  fixed  frequency  will  result  in  less  damping. 
If  the  response  at  that  frequency  was  responsible  for  the  instabil¬ 
ity  then  the  solution  at  the  smaller  time  step  may  actually  be  less 
stable. 

Fig.  10  shows  the  out-of-plane  motion  of  the  free  end  of  the 
chain  for  the  algorithms  that  ran  for  the  full  10  s  at  both  At 
=  0.01  s  and  Af= 0.001  s.  At  Af  =  0.01  s  there  is  little  consistency 
between  the  levels  of  out-of-plane  motion  predicted  by  the  differ¬ 
ent  algorithms.  For  the  solutions  at  Ar=  0.001  s  the  results  for 
out-of-plane  response  appear  roughly  equivalent.  A  trace  of  the 
motion  of  the  free  end  in  the  horizontal  plane  for  Ar  =  0.001  s  and 
\*=-0.3  is  shown  in  Fig.  11.  The  roughly  circular  whirling 
motion  revealed  by  the  trace  after  the  three-dimensional  motion  is 
fully  developed  as  expected  for  this  problem  (Nayfeh  and  Mook 
1995). 

Conclusions 

The  stability  properties  of  the  box,  backward  difference,  trapezoi¬ 
dal  rule,  and  generalized-a  temporal  integration  algorithms  were 
studied.  The  box  method  is  popular  in  cable  dynamics  applica¬ 


tions  because  it  is  second-order  accurate  and  easy  to  implement. 
In  the  harmonically  driven  hanging  chain  problem  considered 
above,  however,  its  poor  numerical  stability  made  solutions  diffi¬ 
cult  or  impossible  to  obtain  beyond  a  certain  point  in  time. 

Backward  differences  have  excellent  numerical  dissipation 
(and  thus  very  good  stability),  but  are  only  first-order  accurate. 
For  the  hanging  chain  problem  this  poor  accuracy  leads  to  nu¬ 
merical  solutions  that  compared  poorly  with  both  analytic  and 
experimental  results.  Unlike  in  the  experiment,  the  simulated 
chain  never  crossed  over  itself.  Trapezoidal  rule  solutions  showed 
good  accuracy,  but  because  of  their  weak  numerical  dissipation, 
relatively  poor  stability  in  the  forced  response  problem.  Of  the 
three  algorithms  that  have  been  popularly  employed  in  cable  dy¬ 
namics  problems  (box,  backward  differences,  trapezoidal  rule)  the 
trapezoidal  rule  appears  to  be  the  best  choice. 

Of  all  the  algorithms  considered,  the  generalized-a  algorithm 
had  the  best  combination  of  accuracy  and  stability.  For  the  har¬ 
monically  driven  chain  it  was  the  only  algorithm  that  produced  a 
simulation  of  the  three-dimensional  whirling  motion  that  develops 
after  crossover  occurs.  While  it  is  slightly  more  complicated  to 
program,  it  has  the  advantage  that  once  it  is  implemented,  box, 
backward  difference,  and  trapezoidal  rule  solutions  can  easily  be 
obtained  through  proper  choice  of  the  parameter  values.  The  best 
choice  for  X30  is  problem  dependent,  but  —  0.5^  X30^  —  0.2  ap¬ 
pears  to  be  a  useful  range  for  many  problems.  Care  must  still  be 
taken  to  insure  adequate  spatial  and  temporal  discretizations  so 
that  the  important  frequency  content  of  the  solution  is  preserved. 
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